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Low-dimensional behavior of large systems of globally coupled oscillators has been intensively investigated 
since the introduction of the Ott-Antonsen ansatz. In this report, we generalize the Ott-Antonsen ansatz to 
second-order Kuramoto models in complex networks. With an additional inertia term, we find a 
low-dimensional behavior similar to the first-order Kuramoto model, derive a self-consistent equation and 
seek the time- dependent derivation of the order parameter. Numerical simulations are also conducted to 
verify our analytical results. 

Synchronization phenomena in large ensembles of coupled systems play a prominent role in many branches 
of natural and social sciences as well as in engineering^'^. The study of collective synchronization has many 
applications including the modeling of the flashing of groups of fireflies^, the collective oscillations of 
pancreatic beta cells^, the human cardiorespiratory system^ and the pedestrian induced oscillations in bridges^. 
A fundamental contribution to the mathematical aspects of collective synchronization was given by Kuramoto^. In 
1975 Kuramoto proposed a model to describe the behaviour of a population of coupled non-linear oscillators, 
employing three key simplifying assumptions^, i.e., (i) the coupling strength was chosen to be homogeneous for all 
pairs of coupled oscillators; (ii) the coupling strength and the natural frequency become finite; and (iii) the number 
of oscillators was considered to be infinite. Diversity in the oscillators properties is usually incorporated by taking 
natural frequencies from a given probability distribution function. The phase transition to synchronization occurs 
when the coupling strength exceeds a threshold, which depends on this probability density function. 

In 2008, Ott and Antonsen^ introduced an ansatz for studying the behaviour of globally coupled oscillators. The 
Ott-Antonsen ansatz has been considered to investigate continuously time -dependent collective behavior^ and 
for the study of delay heterogeneity^^. In addition, such ansatz has enabled to find nonuniversal transitions to 
synchrony in the model with a phase lag for certain unimodal frequency distributions ^\ 

Although these works have provided important contributions to synchronization theory, only oscillators with 
global coupling have been taken into account^"^^. Thus, a natural extension of these works can investigate how 
these results change when different coupling schemes are introduced. Barlev et al.^^ studied the dynamics of 
coupled phase oscillators, but such approach involved integrating AT ordinary differential equations. To overcome 
this limitation, in this report we generalize the Ott-Antonsen ansatz to complex networks in the continuum limit 
to investigate a time- dependent phase transition to synchronization. We reduce the dimension of the system of 
equations from N to the number of possible degrees in the network. 

Motivated by the results of the first- order Kuramoto model, we substantially extend the theory to the second- 
order Kuramoto model. The Kuramoto model with inertia has been widely used for deepening the understanding 
of power grids^^"^^, superconducting Josephson Junctions^^ and many other applications^^'^^. Therefore a theory 
that investigates the low-dimensional character of such systems giving access to their time- dependent behavior 
can bring important new insights into the study of the second-order Kuramoto model. We substantially address 
this problem for what is perhaps the simplest choice of inertia term. In this case, the Fourier series expansion, the 
key approach of the Ott-Antonsen ansatz, no longer applies directly. Thus, a generalized framework for the 
second derivative needs to be developed, as already pointed out in recent studies^^'^°. In order to fill this gap, we 
derive self-consistent equations and seek the time evolution of the order parameter. Comparison of analytical and 
simulation results shows a good agreement. Our results shed light on the impact of the topology on the global 
dynamics. 
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Results 

We consider the first-order Kuramoto model on an unweighted and 
undirected complex network. The state of oscillator / is denoted by its 
phase ^/ ( / = 1 , 2, • • • , AT) , and the governing equation of the modeF is 



dt 



where Q/ stands for the natural frequency of oscillator /, which is 
distributed according to some probability density ^(Q), K specifies 
the homogeneous coupling strength between interconnected nodes, 
and Ay is the element of the adjacency matrix A, i.e., = 1 if nodes / 
and; are connected or = 0, otherwise. 

In uncorrelated networks, if AT approaches infinity (in thermodyn- 
amic limit), the probability of selecting an edge connected to a node 
with degree /c, natural frequency Q, and phase 0 at time t is kP{k)p{k; 
Q, 0, t), where we define P(k) as the degree distribution and p(k; Q, 0, 
t)lk as the probability distribution function of nodes with degree k 
that have natural frequency Q and phase 0 at time t^^^^^-^^ and k the 
average degree. 

To characterize the macroscopic behavior of the oscillators, in the 
continuum limit, we consider the order parameter (see Methods for 
details) 



= dk 



: I dQ^ dOP{k)kp{k; Q, 0, ty^ j | dkP(k)k 



dkP(k)kry'^ 



(2) 



dkP{k)K 



where quantifies the local synchrony of oscillators with degree k 

ry^^^ ^ I I ^^^(^5 0, t)e'^. (3) 

For simplicity, we assume that the natural frequencies Q/ are distrib- 
uted according to an unimodal and symmetric Cauchy-Lorentz dis- 
tribution (g{Q)) (see Methods for details) with zero mean. We set i/^ 
= ij/k = 0 without loss of generality^^. The coupling term in Eq. (1) 

can be written as ^^^^ A^ sin(^^-^,) =kirlm[e'^-^f''''''\ Thus 

the governing Eq. (1) can be rewritten as 



dt 



2i 



(4) 



which shows that the oscillators are coupled via the mean-field order 
parameter r. The restoring force tends to bring each oscillator 
towards equilibrium and the amount of forcing is proportional to 
its degree k. 

The evolution of p(k; Q, 0, t) is governed by the continuity equa- 
tion, i.e., ^ + ^ = 0, where v(k; Q,9,t) = ^. We use the Ott- 

dt dO ^ ' ^ dt 

Antonsen ansatz^ and expand the density function in a Fourier series, 

i.e.. 



p{k; Q, 9, t): 



2k 



1 + 



£[a(^;Q, 



e +C.C. 



(5) 



where c.c stands for the complex conjugate. Substituting the expan- 
sion into Eqs. (3) and in the continuity equation, we get that r^^ = a{k) 
and Yk evolve according to 



drk 
dt 



= -rk+- 



Kkr 



l-rl) for /CG [kmin.krr 



(6) 



where /Cmm and /Cmax are the minimum and the maximum degree, 
respectively. This method works efficiently compared to the ref- 
erence 24 especially when the power law behavior has some cutofP^ 
a{k) therein allows a clear physical interpretation as measuring the 
internal synchrony of the nodes with the same degree k. The global 
order parameter r is a sum of different multiplied by their degree 
and degree distribution (see Eq. (2)). 

To verify the accuracy of the time evolution of the order parameter 
(see Eq. (6)), we compare the time evolution of the order parameter 
r with numerical simulations. Fig. 1 shows the results. Initially, the 
values of oscillators are selected at random from 7i to — 7i, which 
implies that the initial value of each rk{0) tends to zero. In our 
simulations, we set r^iO) = 0.001. As we can see in Fig. 1, the results 
obtained through the solution of the reduced system in Eq. 6 are in 
good agreement with the numerical simulations. 

The analysis above shows the remarkable usefulness of the Ott- 
Antonsen ansatz of the first-order Kuramoto model in complex net- 




Simulation r 
— Analytic r 



1.5 



t 
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Figure 1 | The order parameter as a function of time. Numerical simulations of the Kuramoto model are conducted on a scale-free network (see Methods 
for details). The coupling strength K = 2.5 and 6 are randomly selected from — 7i to 7i at t = 0. 
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Figure 2 | Phases ^and frequencies 6 vs natural frequencies fl, which shows that phase-locked oscillators only exist in red area but not in the 
yellow area. The read area indicates parameter combination of stable fixed point. Stable fixed points and limit cycles coexist in the yellow area. The white 
area represents the existence of limit cycles. The stationary value of the order parameter r could be calculated by simulations or Eq. (11). Thus nodes with 
natural frequencies between ^ — 4\/Kr ^ k, /^j — [ — 3.57, 3.57] are synchronized. The boundary of bistable region are specified by IQI within 
\^4VKr j 7i,Krj = [3.57, 7.18]. 



works, but what happens when we consider the Kuramoto model 
with inertia? The simplest and most straightforward way is to include 
one unity inertia term. This leads to the mean-field character of the 
second-order Kuramoto modeP^"^^ 



~dt 



+ Q/ + iCA;rsin(-^/), 



(7) 



where k varies from the minimal to the maximal degree. 

As shown in Eq. (5), the main idea of the Ott-Antonsen ansatz is to 
expand the probability density p{k; Q, 0, t) in a Fourier series in 0. 
For the Kuramoto model with inertia, the probability density 

is also a function of the additional term 0. As 0 varies 

from —00 to it is not possible to follow the same procedure to 
derive the nonlinear evolution of the order parameter r. Due to the 
existence of the inertia term and the bistable area of the stability 
diagram^^, we rewrite Eq. (6) with two functions A{K k) and/(iC /c, 
r) = a(KA;)Vand get 



Melnikov's method^^ is used to show that oscillators are within a 
stable fixed point area as j5 ^ 0 and 1 < 4j5/7i; only limit-cycle 
oscillators exist for J > 1; limit cycles and stable fixed points coexist 
otherwise. 

Let us first investigate the stationary states of phases 0 and 0 in 
terms of the natural frequencies Q separately. In Fig. 2, every single 
point represents the state of one oscillator at time T(T>1) using 
simulations with N = 10000 nodes and degree K = 10. It is interest- 
ing to find that instead of three different regions mentioned above, 
the oscillators fall into either of the following two groups, (i) If the 
natural frequencies of nodes are within the boundary of the phase 

synchronization regime [^^lower,^^upper] = ^ — 4:\/Kr ^ tt, aVKt ^ Trj 

which is the same as the above stable fixed points area, these nodes 
converge to fixed points and the stationary state of phases are func- 
tions of Q, which are equal to arcsin(Q/(iCr)). This boundary is 
smaller than that of the Kuramoto model, in which oscillators are 
in the locked state for all IQI < Kr^^. (ii) In contrast, the oscillators 



dvk 
dt 



-n+A{Kk)- 



r 1- 



(8) 



where A{K k) indicates the effective coupling strength and a, h and c 
are constant. /(iC /c, r) is a high-order term and is used to adjust the 
stationary solution. For the Kuramoto model without inertia, we get 
A{Kk) = Kk eiud a = 0. 

In order to solve Eq. (8), we first investigate the nonlinear 
dynamics on fully connected networks. In this case, we normalize 
the coupling strength from KNtoK. For the sake of convenience, we 
change the time scale to t = VKrt, which yields 



d% 

dT^ 



^dOi^ 
di 



+ J/ + sin(-^/), 



(9) 



where P=\ j yKr and J/ = Qi/(Kr). Thus P is identical for all oscil- 
lators and the diversity of If is due to its natural frequency. According 
to the parameter space^^'^^, nodes are divided into three groups. 



with |Q| >\\/Kr j n are drift. Thus, in networks, instead of three 

different areas of single pendulum model, only two distinct areas 
could exist: fixed point and limit cycle. Nodes with the same natural 
frequency are either converging to single fixed points or oscillating 
periodically; and nodes always return to previous states even after 
large perturbations. 

To investigate how the phase synchronization boundary changes 
with different coupling strengths, we project the Fig. 2 on the I-p 
parameter space and color the oscillators according to their station- 
ary states in the parameter space. A comparison between the 
dynamics with average degree 10 and that with 30 is shown in 
Fig. 3. We can see that oscillators with the same coupling share the 
same P axis and the diversity of I is due to the distribution of the 
natural frequencies Q. All synchronized nodes are inside the syn- 
chronized area, which is at the right side of the line I = 4p/n. 

Therefore, after substituting the boundaries of the synchronized 
natural frequencies [Qiower^ ^upper] and the Cauchy-Lorentz distri- 
bution into the definition of the order parameter r. 
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Figure 3 | The definitions of three shaded areas are the same as that in Figure 2. Two boundaries are compared between coupling strengths 10 and 30. If 
oscillators are in locked state with black color and with Chartreuse color otherwise. Increasing the coupling strength iC further, the vertical line moves to 
the left. 



Jq,„ 



(10) 



model (proportional to 4\/k/ tt), we set A{K)=4\/k / n. When 



r = 0, 



where Os denotes the synchronized oscillator sin (0^) = I. Performing 
some mathematical manipulations, we get 



/(iC, r) = r-4VKr{l-r^) j (In), 



(12) 



2 

71^ 



A^Kr+{Krf 
^(^ 



-\6Kr 



(11) 



Due to the difference of boundaries between the first- order Kura- 



and this stationary solution should be met by the self- consistent Eq. 
(11). Here, we use numerical methods to calculate the values of a, h 
and c. As shown in Fig. 4, after substituting the stationary solutions K 
and r of Eq. (1 1) into Eq. (12),/(K, r) is colored in red and we get the 
values a = 0.389, b= 1/4 and c= 3. When r is small,/(K, r) is close to 
0 and cannot influence the time evolution of the order parameter r(t). 



moto model (proportional to K) and the second-order Kuramoto or vary the stationary solution, otherwise. 
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Figure 4 | /(IC, r) as a function of stationary solution of self-consistent equation colored in red and the fitting curve colored in blue. 
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Figure 5 | Order parameter rvs coupling strengths Kin scale-free networks (see Methods for details). The red curves indicate the results from 
simulations on the same network as in Figure 1. For each coupling, initial values of 6 randomly select from [ — 7i, 7i] and we set 6 = 0. The green dots shows 
analytic prediction of the stationary r(f) based on the self-consistent Eq. (8). 



Let us consider again the nonlinear evolution of the order para- 
meter r in complex networks. From the above analysis, we get that 

A{Kk) =4VKk j K. To check the validity of this assumption, we 

compare the stationary solution with simulation results in Fig. 5. 
The theoretical predictions (green lines derived from Eq. (8) with 
effective coupling and f{K k, r)) are in agreement with red lines of 
numerical simulations. 

The nonlinear evolution of r{t) is illustrated in Fig. 6, for a selection 
of coupling strengths K. Initial values of Oi and Oi are the same as in 
Fig. 5. For the order parameter formulation the initial value of r is set 
to a small value (r(0) 1). The r formulation of Eq. (8) does not 



only reproduce the stationary states in Fig. 5, but also matches the 
transition to synchrony. The analytic results and simulation results 
are in good agreement. 

Conclusions 

In conclusion, we proposed a generalization for the Ott-Antonsen 
ansatz to complex networks with a Cauchy-Lorentz distribution of 
the natural frequency for the Kuramoto model. Compared to the 
ensemble approach^^, the dimension of ordinary differential equa- 
tions was reduced from N to the number of possible degrees in the 
network. We have investigated the collective dynamics of the 
Kuramoto model with inertia and found the synchronization bound- 



1.0 




t 



Figure 6 | Order parameter r{t) vs time tin scale-free networks (see Methods for details). The simulations are conducted on the same network 

and the coupling strength K = \ and iC = 3. Blue and yellow dots are analytic results got from Eq. (8). In simulations, initial values of 0 are randomly 

selected from —%to% and that of 0 close to 0. 
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ary is [ — 4:\/Kr j n, 4:\/Kr j n] instead of [—Kr, Kr] as in the 

Kurmoto model without inertia. Based on these results, we analyt- 
ically derived self-consistent equations for the order parameter and 
nonlinear time-dependent order parameter. The agreement between 
the analytical and simulation results is excellent. 

Methods 

The networks. The model has been implemented on undirected and unweighted 
scale-free networks with N = 10000, P{k) oc and /c > 5. 

Numerical integration. Eqs. (6) and (12) are solved by a order Runge-Kutta 
method with time step h = 0.01 and with the Cauchy-Lorentz distribution 

Order parameter. In complex networks, in order to understand the dynamics of the 



system, it is natural to use the definition of order parameter r^^ as re^^ 



instead of the definition r 



N 



- , which accounts for the mean-field in the fully 



connected graph regime. 

The magnitude r ^ [0, 1] quantifies the phase coherence, while \J/ denotes the 
average phase of the system. In particular, r ~ 0, if the phases are randomly distributed 
over [0, In] and all nodes oscillate at its natural frequency. On the other hand, if all 
oscillators run as a giant component, r~ 1. The system is known to exhibit a phase 
transition from the asynchronous state (r^O) to the synchronous one (r~ 1) at a 
certain critical value /l^ characterizing the onset of partial synchronization and, for 
unimodal and symmetric frequency distributions g{Q), the transition is continuous. It 

2 

turns out that for uncorrelated networks, /l^ is given by /Ic = ^ where /Imax 



is the maximal eigenvalue of the adjacency matrix. 
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